Parameterized report
Package installation
This notebook is a template you can use for evaluating your forecaster against submissions from COVIDhub. While the template is for forecasting for state deaths, it can be modified easily to forecast for county cases. We will point out what needs to be modified for counties in the template (other than the obvious variable names and R chunk names).
First, specify the dates we want forecasts for. n_dates - 4 is the most recent date with ground truth for 4 weeks ahead, so that is usually the last date you want to get forecasts for. To cover a longer time period, we usually have the forecast dates spaced 2-4 weeks apart.
our_pred_dates <- get_covidhub_forecast_dates("CMU-TimeSeries")
n_dates <- length(our_pred_dates)
# n_dates - 4 is often the most recent date with ground truth for 4 weeks ahead
# dates spaced out 2 weeks ahead to see different behaviors
forecast_dates <- our_pred_dates[n_dates- 2 *5:2]To use this template for your forecaster, Your main job is to specify the arguments to evalcast::get_predictions, which are listed below. See ?evalcast::get_predictions for details on what these arguments should be; I give some further notes below.
forecastername_of_forecastersignalsforecast_dates: Just useforecast_datesas defined in the chunk above.incidence_period: In our setting, almost always"epiweek".apply_corrections: These corrections are at the daily level. If you want to use this argument, I suggest extracting the relevant function/code from the production parameter script (utilities/production_params.R).forecaster_args
The code below is an example of what these arguments might be.
aheads <- as.numeric(params$weeks)
lags <- c(0, 7, 14)
state_forecaster <- larrys_anteater
state_forecaster_name <- "larry's anteater"
# don't put an as_of column in the signals dataframe: this will be
# taken care of correctly by evalcast::get_predictions
# (i.e. as_of = the relevant forecast date)
state_forecaster_signals <- dplyr::tibble(
data_source = "jhu-csse",
signal = c("deaths_incidence_num",
"confirmed_incidence_num"
),
start_day = lubridate::ymd("2020-03-01"),
geo_type = "state"
)
state_forecaster_args <- list(
ahead = aheads,
lags = lags,
tau = evalcast::covidhub_probs(), # 23 quantiles
training_window_size = 90,
featurize = animalia::make_7dav_featurizer()
)State correction parameters (taken directly from production_params.R):
state_corrections_params <- zookeeper::default_state_params(
# many other options, see the function documentation
data_source = state_forecaster_signals$data_source,
signal = state_forecaster_signals$signal,
geo_type = state_forecaster_signals$geo_type)
state_corrector <- zookeeper::make_state_corrector(
params = state_corrections_params,
# data, locations, times to do special correction processing
manual_flags = tibble::tibble(
data_source = "jhu-csse",
signal = c(rep("deaths_incidence_num", 3),
"confirmed_incidence_num",
## from JHU-CSSE notes 2021-04-17, 2021-04-18
"deaths_incidence_num",
"deaths_incidence_num",
"confirmed_incidence_num",
"confirmed_incidence_num",
## ## from JHU-CSSE notes 2021-04-24, 2021-04-25
## "confirmed_incidence_num"
## from JHU-CSSE notes 2021-05-02
"deaths_incidence_num"
),
geo_value = c("va","ky","ok","ok",
## from JHU-CSSE notes 2021-04-17, 2021-04-18
"ak","mi","mo","al",
## ## from JHU-CSSE notes 2021-04-24, 2021-04-25
## "al"
## from JHU-CSSE notes 2021-05-02
"wv"
),
time_value = list(
seq(lubridate::ymd("2021-02-21"), lubridate::ymd("2021-03-04"), by = 1),
lubridate::ymd(c("2021-03-18","2021-03-19")),
lubridate::ymd("2021-04-07"),
lubridate::ymd("2021-04-07"),
## from JHU-CSSE notes 2021-04-17, 2021-04-18
lubridate::ymd("2021-04-15"),
lubridate::ymd(c("2021-04-01", "2021-04-03", "2021-04-06", "2021-04-08", "2021-04-10", "2021-04-13", "2021-04-15", "2021-04-17",
## seems to be ongoing as of week of 2021-04-27
"2021-04-20", "2021-04-22", "2021-04-24",
"2021-04-27", "2021-04-29", "2021-05-01"
)),
lubridate::ymd("2021-04-17"),
lubridate::ymd("2021-04-13","2021-04-20"),
## ## from JHU-CSSE notes 2021-04-24, 2021-04-25
## lubridate::ymd("2021-04-20")
## from JHU-CSSE notes 2021-05-02
lubridate::ymd("2021-04-27")
),
max_lag = c(rep(90, 4),
## from JHU-CSSE notes 2021-04-17, 2021-04-18
75, 150, 150, 180,
## from JHU-CSSE notes 2021-04-24, 2021-04-25
## as.integer(as.Date("2021-04-20") - as.Date("2020-10-23"))
## from JHU-CSSE notes 2021-05-02; just assign an arbitrary large value due to lack of accessible details
180
)
)
)Get the predictions for our forecaster (you can leave this call as-is, unless you want to change incidence_period):
set.seed(111) # for reproducibility as corrections have randomness
state_predictions <-
evalcast::get_predictions(
forecaster = state_forecaster,
name_of_forecaster = state_forecaster_name,
signals = state_forecaster_signals,
forecast_dates = forecast_dates,
incidence_period = "epiweek",
apply_corrections = state_corrector,
forecaster_args = state_forecaster_args)Sometimes you might want to post-process the predictions: you can do them here. In this example, our model could give us negative predictions and we choose to enforce a non-negativity constraint.
state_predictions$value <- pmax(state_predictions$value, 0)Get the competition and evaluate the predictions on the error metrics:
competition <- params$forecasters
if ("Karlen-pypm" %in% competition) {
competition <- competition[-which(competition == "Karlen-pypm")]
}
if (!"COVIDhub-baseline" %in% competition) {
competition[[length(competition)+1]] <- "COVIDhub-baseline"
}
submitted <- lapply(competition[1:length(params$forecasters)-1], get_covidhub_predictions,
forecast_dates = forecast_dates,
signal = "deaths_incidence_num",
ahead = aheads)
if ("Karlen-pypm" %in% params$forecasters) {
submitted[[length(params$forecasters)]] <- get_covidhub_predictions("Karlen-pypm",
forecast_dates = forecast_dates - 1,
signal = "deaths_incidence_num") %>%
mutate(forecast_date = forecast_date + 1)
}submitted <- bind_rows(submitted) %>% filter(ahead < 5)
all_predictions <- bind_rows(state_predictions, submitted)# for county prediction, set geo_type = "county"
results <- evaluate_covid_predictions(all_predictions,
backfill_buffer = 0,
geo_type = "state") %>%
intersect_averagers(c("forecaster"), c("forecast_date", "geo_value"))
results %>%
group_by(forecast_date) %>%
summarise(n_distinct(geo_value))## # A tibble: 3 x 2
## forecast_date `n_distinct(geo_value)`
## <date> <int>
## 1 2021-05-24 48
## 2 2021-06-07 48
## 3 2021-06-21 48##without using API
# predictions_cards <- readRDS("~/Documents/DSSG 2021/dssg_covidcast/predictions_cards.rds")
# results <- evaluate_covid_predictions(predictions_cards,
# backfill_buffer = 0,
# geo_type = "state") %>%
# intersect_averagers(c("forecaster"), c("forecast_date", "geo_value"))Overall AE, WIS, Coverage 80
NOTE: Results are based on the following numbers of common locations
weeks.label = c("1 week ahead", "2 weeks ahead", "3 weeks ahead", "4 weeks ahead")
names(weeks.label) = c(1, 2, 3, 4)
subtitle = sprintf("Forecasts made over %s to %s",
format(min(forecast_dates), "%B %d, %Y"),
format(max(forecast_dates), "%B %d, %Y"))
plot_ae <-
plot_canonical(results,
x = "ahead",
y = "ae",
aggr = mean) +
labs(title = subtitle, x= "Weeks Ahead", y = "Mean AE",color='Forecasters') +
geom_point(aes(text=sprintf("Weeks Ahead: %s<br>Average Error: %s <br>Forecaster: %s",
ahead,
round(ae, digits=2),
color)),
alpha = 0.05) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5)) +
scale_y_log10()
if (params$colorblind_palette) {
plot_ae <- plot_ae +
scale_color_viridis_d()
}
ggplotly(plot_ae, tooltip="text", width=1000) %>%
layout(hoverlabel = list(bgcolor = "white"))plot_wis <-
plot_canonical(results,
x = "forecast_date",
y = "wis",
aggr = mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead") +
labs(title = subtitle,
x = "Forecast Dates",
y = "Mean WIS",
color = "Forecasters") +
geom_point(aes(text=sprintf("Forecast Date: %s<br>Mean WIS: %s <br>Forecaster: %s",
forecast_date,
round(wis, digits = 2),
color)),
alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead=weeks.label)) +
theme(strip.background = element_rect(fill = "white")) +
theme_minimal() +
theme(plot.title = element_text(hjust = "center")) +
scale_y_log10()
if (params$colorblind_palette) {
plot_wis <- plot_wis +
scale_color_viridis_d()
}
ggplotly(plot_wis, tooltip="text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))plot_cover <-
plot_canonical(results,
x = "forecast_date",
y = "coverage_80",
aggr = mean,
grp_vars = c("forecaster","ahead"),
facet_rows = "ahead") +
labs(title = subtitle, x= "Forecast date", y = "Mean Coverage 80", color='Forecasters') +
geom_point(aes(text = sprintf("Forecast Date: %s<br>Coverage: %s <br>Forecaster: %s",
forecast_date,
round(coverage_80, digits = 2),
color)),
alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = weeks.label)) +
theme(strip.background = element_rect(fill = "white")) +
theme_minimal() +
theme(plot.title = element_text(hjust = "center")) +
geom_line(mapping = aes(y = .8))
if (params$colorblind_palette) {
plot_cover <- plot_cover +
scale_color_viridis_d()
}
ggplotly(plot_cover, tooltip="text", height=800, width=1000) %>%
layout(hoverlabel = list(bgcolor = "white"))(Geometric) Mean relative WIS
Relative to baseline; scale first then take the geometric mean, ignoring a few 0’s. I think this is potentially more useful than the median/mean for relative WIS (or relative AE), but I haven’t completely thought it through. Putting the results here to be provocative.
geom_mean <- function(x) prod(x)^(1/length(x))
#geom_mean <- exp(mean(log((x+1)/(y+1)))) #still need to figure this out
mean_wis <-
plot_canonical(results %>%
filter(wis > 0),
x = "ahead",
y = "wis",
aggr = geom_mean,
base_forecaster = "COVIDhub-baseline",
scale_before_aggr = TRUE) +
labs(title = subtitle,
x = "Weeks Ahead",
y = "Mean (geometric) relative WIS",
color = "Forecasters") +
geom_point(aes(text = sprintf("Weeks Ahead: %s<br>WIS: %s <br>Forecaster: %s",
ahead,
round(wis, digits = 2),
color)),
alpha = 0.05) +
theme_minimal() +
theme(plot.title = element_text(hjust = "center")) +
geom_line(mapping = aes(y = 1))
if (params$colorblind_palette) {
mean_wis <- mean_wis +
scale_color_viridis_d()
}
ggplotly(mean_wis, tooltip="text", width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))mean_wis_forecast_date <-
plot_canonical(results %>%
filter(wis > 0),
x = "forecast_date",
y = "wis",
aggr = geom_mean,
facet_rows = "ahead",
grp_vars = c("forecaster", "ahead"),
base_forecaster = "COVIDhub-baseline",
scale_before_aggr = TRUE) +
theme(legend.position = "bottom") +
labs(title = subtitle,
x = "Forecast date",
y = "Mean (geometric) relative WIS",
color = "Forecasters") +
geom_point(aes(text = sprintf("Forecast Date: %s<br>WIS: %s <br>Forecaster: %s",
forecast_date,
round(wis, digits = 2),
color)),
alpha = 0.05) +
facet_wrap(~ahead, nrow = 4, labeller = labeller(ahead = weeks.label)) +
theme(strip.background = element_rect(fill = "white")) +
theme_minimal() +
theme(plot.title = element_text(hjust = "center")) +
geom_line(mapping = aes(y = 1))
if (params$colorblind_palette) {
mean_wis_forecast_date <- mean_wis_forecast_date +
scale_color_viridis_d()
}
ggplotly(mean_wis_forecast_date, tooltip = "text", height=800, width= 1000) %>%
layout(hoverlabel = list(bgcolor = "white"))Scores by target date (not forecast date)
Maps (mean score over forecast dates and aheads)
Note that the results are scaled by population.
library(sf)
maps <- results %>%
group_by(geo_value, forecaster) %>%
summarise(across(wis:ae, mean)) %>%
left_join(animalia::state_population, by = "geo_value") %>%
mutate(across(wis:ae, ~ .x / population * 1e5)) %>%
pivot_longer(wis:ae, names_to = "score") %>%
group_by(score) %>%
mutate(time_value = Sys.Date(),
r = max(value)) %>%
group_by(forecaster, .add = TRUE) %>%
group_split()
# for county prediction, set geo_type = "county"
maps <- purrr::map(maps,
~as.covidcast_signal(
.x, signal = .x$score[1],
data_source = .x$forecaster[1],
geo_type = "state"))
maps <- purrr::map(maps,
~plot(.x,
choro_col = scales::viridis_pal()(3),
range = c(0,.x$r[1])))
nfcasts <- length(unique(results$forecaster))Mean AE
# original code
cowplot::plot_grid(plotlist = maps[1:nfcasts], ncol = 3)
Mean WIS
cowplot::plot_grid(plotlist = maps[(nfcasts+1):length(maps)], ncol = 3)
Trajectory plots
For county predictions, you might want to change the fig.height and fig.width chunk options here (in other notebooks we use fig.height = 120 and fig.width = 30).
# for county prediction, set geo_type = "county"
state.label = c(state.name, "Washington D.C.", "Porto Rico")
names(state.label) = c(tolower(state.abb), "dc", "pr")
pd <- evalcast:::setup_plot_trajectory(
bind_rows(state_predictions, submitted %>%
filter(forecaster == "CMU-TimeSeries")),
intervals = 0.8,
geo_type = "state",
start_day = min(forecast_dates) - 60)
g <- ggplot(pd$truth_df, mapping = aes(x = target_end_date))
# build the fan
g <- g + geom_ribbon(
data = pd$quantiles_df %>%
mutate(upper = round(upper, 2),
lower = round(lower, 2)),
mapping = aes(ymin = lower,
ymax = upper,
fill = forecaster,
group = interaction(forecaster, forecast_date)),
alpha = .1) +
scale_fill_viridis_d(begin=.15, end=.85)
# line layer
g <- g +
#geom_line(aes(y = .data$value.y), color = "#3182BD") + # corrected
geom_line(aes(y = value), size = .5) + # reported
geom_line(data = pd$points_df,
mapping = aes(y = value,
color = forecaster,
group = interaction(forecaster, forecast_date)),
size = .5) +
geom_point(aes(y = value), size = 1) + # reported gets dots
geom_point(data = pd$points_df %>%
mutate(value = round(value, 2)),
mapping = aes(y = value, color = forecaster),
size = 1) +
scale_color_viridis_d(begin=.15, end=.85)states <- g +
facet_wrap(~geo_value,
scales = "free_y",
ncol = 3,
labeller = labeller(geo_value = state.label)) +
labs(x = "", y = "") +
theme_bw() +
theme(legend.position = "none", strip.background = element_rect(fill = "white"))
ggplotly(states, height=5000, width= 1000)